

havesave = 1b
undefine,rb3
native = 0b

IF ~ havesave THEN fill_in_rb3_grs,/CLEM


;; Get galactic params
defsysv, '!MW', exists = exists
IF NOT exists THEN galactic_params 

restore,'./irdc_dist_model/bgps_rb3.sav',/ver

IF native THEN ind = where(rb3.grs EQ 1, nind) ELSE BEGIN
   restore,'./irdc_dist_model/new_grs_look.sav',/ver
   ind = where_array(new_grs,rb3.cnum,nind)
ENDELSE

cnum = rb3[ind].cnum
rb3 = rb3[ind]
print,nind

root = './local/output/grs_hi'
post = './local/postage/'
ubc  = './local/ubc_glimpse_proc/smooth30arc'

gamma = FLUX2TAU_BGPS(20.d,/mjy)

;; Also, write out text file for handwriting...
openw,lun,'./local/output/GRS_HI_Verification.txt',/GET_LUN

hisa_int = fltarr(nind)

FOR i=0L, nind-1 DO BEGIN
   
   scnum = string(cnum[i],format="(I04)")
   grs = readfits(post+'grs'+scnum+'.fits',ghdr,/silent)
   hi  = readfits(post+'hisa'+scnum+'.fits',hhdr,/silent)
   lab = readfits(post+'label'+scnum+'.fits',lhdr,/silent)
   smoo = readfits(ubc+scnum+'.fits',shdr,/silent)
   bgps_smoo = readfits(post+'smooth'+scnum+'.fits',bshdr,/silent)
   
   plot_wcs_axes,lhdr,10,lonarr=xc,latarr=yc 
   
   ;; Weight Mask -- derived from BGPS Label image
   undefine,wtmask
   wtmask  = aperture_wtmask(lhdr, 40., rb3[i].l, rb3[i].b)
      
   IF total(size(grs,/DIM)) EQ 0 OR total(size(hi,/DIM)) EQ 0 THEN BEGIN
      print,'NO POSTAGE!  ',rb3[i].cnum, rb3[i].l, rb3[i].b, rb3[i].vlsr
   ENDIF ELSE BEGIN
      myps,root+scnum+'.eps',xsize=18,ysize=6
      
      multiplot,[3,1],xgap=0.04,mtitle='BGPS #'+scnum
      
      ;;============================
      ;; GRS Panel
      undefine,grsr,hisar,smoor
      plotmap,grs,ghdr,ct=0,title='GRS',range=grsr
      cgLoadct,40,/silent
      contour,bgps_smoo,xc,yc,levels=[.025,.05,0.1,0.25,0.5],$
              c_colors=[1,2,3,4,5]*50,/overplot,thick=1.5
      cgLoadct,0,/silent
      contour,lab,xc,yc,level=[0.5],color=cgColor('Red'),thick=4,/overplot
      contour,wtmask,xc,yc,levels=[0.5],thick=5,color=cgColor('Dodger Blue'),$
              /over,c_linestyle=2
      plotsym,3,2.0,/fill
      plots,rb3[i].grs_l,rb3[i].grs_b,psym=8,color=cgColor('Bisque')
      plotsym,3,2.0
      plots,rb3[i].grs_l,rb3[i].grs_b,psym=8,color=cgColor('Opposite')
      ;; plots,rb3[i].grs_l,rb3[i].grs_b,psym=2,thick=4,color=cgColor('TAN7'),$
      ;;       symsize=2
      
      multiplot
      
      ;;============================
      ;; HISA Panel
      plotmap,hi,hhdr,ct=0,title='HISA',range=hisar
      wtmask2     = aperture_wtmask(hhdr, 40., rb3[i].l, rb3[i].b)
      hisa_int[i] = total(wtmask2*hi,/NAN) / total(wtmask2)
      cgLoadct,40,/silent
      contour,bgps_smoo,xc,yc,levels=[.025,.05,0.1,0.25,0.5],$
              c_colors=[1,2,3,4,5]*50,/overplot,thick=1.5
      cgLoadct,0,/silent
      contour,lab,xc,yc,level=[0.5],color=cgColor('Red'),thick=4,/overplot
      contour,wtmask,xc,yc,levels=[0.5],thick=5,color=cgColor('Dodger Blue'),$
              /over,c_linestyle=2
      plotsym,3,2.0,/fill
      plots,rb3[i].grs_l,rb3[i].grs_b,psym=8,color=cgColor('Bisque')
      plotsym,3,2.0
      plots,rb3[i].grs_l,rb3[i].grs_b,psym=8,color=cgColor('Opposite')
      ;; plots,rb3[i].grs_l,rb3[i].grs_b,psym=2,thick=4,color=cgColor('TAN7'),$
      ;;       symsize=2
      
      multiplot
      
      ;;============================
      ;; GLIMPSE Panel
      plotmap,smoo,shdr,ct=0,title='Smoothed GLIMPSE',range=smoor
      cgLoadct,40,/silent
      contour,bgps_smoo,xc,yc,levels=[.025,.05,0.1,0.25,0.5],$
              c_colors=[1,2,3,4,5]*50,/overplot,thick=1.5
      cgLoadct,0,/silent
      contour,lab,xc,yc,level=[0.5],color=cgColor('Red'),thick=4,/overplot
      contour,wtmask,xc,yc,levels=[0.5],thick=5,color=cgColor('Dodger Blue'),$
              /over,c_linestyle=2
      plotsym,3,2.0,/fill
      plots,rb3[i].grs_l,rb3[i].grs_b,psym=8,color=cgColor('Bisque')
      plotsym,3,2.0
      plots,rb3[i].grs_l,rb3[i].grs_b,psym=8,color=cgColor('Opposite')
      ;; plots,rb3[i].grs_l,rb3[i].grs_b,psym=2,thick=4,color=cgColor('TAN7'),$
      ;;       symsize=2
      
      multiplot,/reset
      multiplot,/default
      multiplot,/reset
      
      cgLoadct,0,/silent
      cgColorbar,/vert,range=grsr,position=[0.05,0.10,0.06,0.90],/yst,$
                 title='GRS Integrated Intensity [K km/s]',charsize=0.7,$
                 format="(F0.2)"
      
      cgColorbar,/top,range=hisar,position=[0.56,0.91,0.86,0.93],/xst,$
                 title='HISA Integrated Intensity (On-Off) [K km/s]',$
                 charsize=0.7,format="(F0.2)",/right
      
      cgColorbar,/vert,range=smoor,position=[0.95,0.10,0.96,0.90],/yst,$
                 title='GLIMPSE Intensity [MJy sr!u-1!n]',$
                 charsize=0.7,format="(F0.2)",/right
      
      
      tand = !MW.R0*cos(rb3[i].l * !dtor)/1.d3
      fdata = 1.d - (rb3[i].c_meas / (1 - exp(-gamma*rb3[i].s_peak)))
      
      
      str = 'd!dGRS!n = '+string(rb3[i].grs_dist,format="(F0.2)")+' kpc    '+$
            'd!dtan!n = '+string(tand,format="(F0.2)")+' kpc    '+$
            'd!dmodel!n = '+string(rb3[i].fitdist,format="(F0.2)")+' kpc    '+$
            'd!dmorph!n = '+string(rb3[i].morph_d,format="(F0.2)")+' kpc    '+$
            'V!dLSR,GRS!n = '+string(rb3[i].grs_v,format="(F0.1)")+' km/s    '+$
            'V!dLSR,BGPS!n = '+string(rb3[i].vlsr,format="(F0.1)")+' km/s    '+$
            'C = '+string(rb3[i].c_meas,format="(F0.3)")+'    '+$
            'S!d1.1!n = '+string(rb3[i].s_peak/1.d3,format="(F0.2)")+' Jy   '+$
            'f!ddata,20!n = '+string(fdata,format="(F0.2)")
      cgText,/norm,0.5,0.05,str,align=0.5,charsize=1.0
      
      myps,/done
      
      str = scnum+'  '+string(rb3[i].grs_dist,format="(F0.2)")+'  '+$
            string(tand,format="(F0.2)")+'  '+$
            string(rb3[i].fitdist,format="(F0.2)")+'  '+$
            string(rb3[i].morph_d,format="(F0.2)")
      printf,lun,str
      
   ENDELSE
   
ENDFOR

close,lun
free_lun,lun

myps,'./irdc_dist_model/analysis_plots/hisa_vs_grsdist.eps',xsize=10,ysize=4.5

multiplot,[2,1],xgap=0.04

plot,rb3.grs_dist,hisa_int,psym=8,xtit='d!dGRS!n [kpc]',$
     ytit='I!dHISA!n [K km/s]'

linfit_plot,rb3.grs_dist,hisa_int,linestyle=2,color=cgColor('Red'),thick=2,$
            /legend,linsize=0.6,charsize=0.8

multiplot,/doxaxis,/doyaxis

tand = !MW.R0*cos(rb3.l * !dtor)/1.d3

plot,rb3.grs_dist/tand,hisa_int,psym=8,xtit='d!dGRS!n / d!dtan!n',$
     ytit='I!dHISA!n [K km/s]'

linfit_plot,rb3.grs_dist/tand,hisa_int,linestyle=2,color=cgColor('Red'),$
            thick=2,/legend,linsize=0.6,charsize=0.8


plots,1.1,!y.crange[0]
plots,1.1,-5,/continue,linestyle=3,color=cgColor('YGB5'),thick=2
plots,!x.crange[1],-5,/continue,linestyle=3,color=cgColor('YGB5'),thick=2

myps,/done,/mp


wtfind = WHERE(rb3.grs_dist/tand GE 1.1 AND hisa_int LE -5, nwtf)
print,nwtf
print,rb3[wtfind].cnum

END
